Comparing numerical and analytical calculations of post-ISCO ringdown amplitudes 
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We numerically compute the ringdown amplitudes following the plunge of a particle from the 
innermost stable circular orbit (ISCO) of a Schwarzschild black hole in the extreme-mass ratio 
limit. We show that the ringdown amplitudes computed in this way are in good agreement with a 
recent analytical calculation [l|]. 
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PACS numbers: 04.40.Dg, 04.62. +v, 95.30.Sf 

I. INTRODUCTION 

In this paper we will study a compact object plunging 
into a much more massive, nonrotating (Schwarzschild) 
black hole from the innermost stable circular or- 
bit (ISCO), located (in Schwarzschild coordinates) at 
nsco = 3rs, where rg = 2M is the Schwarzschild radius 
and M is the black hole mass. This post-ISCO plunge 
trajectory is of special interest because it is "universal" . 
It has long been known 0] that the eccentricity of bod- 
ies orbiting a black hole must decrease in the Newtonian 
regime of low velocities and large separations (r ^ Vg) 
during a gravitational- wave driven inspiral . In this sense, 
a plunge from a quasicircular ISCO represents a "Kep- 
lerian attractor" : for long, gravitational- wave driven in- 
spirals the eccentricity should be essentially zero by the 
time the particle reaches the ISCO. However, some as- 
trophysical scenarios do predict the possibility of orbits 
retaining nonzero eccentricity all the way down to plunge 
(see e.g. H^). 

This "universal" plunge trajectory for particles falling 
from a quasicircular ISCO leaves a very specific signature 
in the quasinormal ringing of the final black hole. The 
amplitudes of the quasinormal modes excited in the pro- 
cess were computed in Ref. In this paper we confirm 
the predictions of that paper by computing gravitational 
radiation with a frequency-domain perturbative code de- 
veloped and tested in Ref. Q , and we verify that the two 
calculations are in very good agreement. Therefore any 
model for extreme mass ratio inspirals leading to plunge 
from a quasicircular ISCO should match the ringdown 
signal predicted in Ref. [l| at late times. It will be in- 
teresting to generalize the present results to comparable 
mass ratio binaries. 

The paper is organized as follows. In section |TT] we 
describe the post-ISCO plunge trajectory and the nu- 
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merical algorithm to compute the gravitational radiation 
produced by the plunging particle, along with the result- 
ing waveforms. In section IIIII we extract the ringdown 
amplitudes from these waveforms and present a compar- 
ison with the results of Ref. [l| . 



II. RADIATION SOURCED BY INFALLING 
OBJECT 

The post-ISCO plunge trajectory is described by the 
coordinates t, r,(j) as a function of the proper time ^ 
(without loss of generality we take the orbit to lie i 
the equatorial plane, i.e. = 7r/2). The trajectory 
solution of the geodesic equations 
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^ = Tr 
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where /(r) = 1 — Tg/r and the energy and angular mo- 
mentum per unit mass have the values corresponding to 
a particle at the ISCO: 
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The plunge trajectory can be written down analytically: 
cf. Eq. (2.7)-(2.11) and Fig. 1 of Ref. In perturbation 
theory, the gravitational radiation at infinity can be de- 
termined from a knowledge of the Sasaki-Nakamura wave 
function Xim [3| (for the extensive literature on gravita- 
tional waves from particles falling into black holes, see 
0, @ and Appendix C of [l^ ) . In the frequency domain, 
the Sasaki-Nakamura equation can be written in the form 
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Here (Z, m) are (tensor) spherical harmonic indices result- 
ing from a separation of the angular variables, uj is the 
Fourier frequency of the perturbation and A = r(r— 2A/). 
The boundary conditions dictate that we should have 
outgoing waves at infinity and ingoing waves at the BH 
horizon: 
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The source term Sim in the Sasaki-Nakamura equation 
(O is determined by the point-particle trajectory, and it 
can be found in Ref . [1] . In terms of the Sasaki-Nakamura 
wavefunctions, the plus- and cross- polarization ampli- 
tudes are given by 
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In our comparisons we will always consider, for simplicity, 
the outgoing amplitudes of the Sasaki-Nakamura wave- 
function in the time domain, i.e. 
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Numerically, the problem is to determine the (com- 
plex) amplitudes Xf^ (t) for plunge trajectories with en- 
ergy and angular momentum given by Eq. Q. In order 
to start the plunge, we must displace the particle from 
the ISCO location by a small quantity e: 



''o = nsco (1 - e) 
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We will present comparisons for e = 5 x 10^""* and e = 
10~^, but we verified that our results arc robust by com- 
puting the radiation for several other values of e, includ- 
ing e = 5 X 10-"*, 10-3, 5 X 10-^, IQ-^, 5 x 10"^, IQ-^. 

Our numerical integrations use a modification of the 
C++ program described in Ref. [1], and we refer the 
reader to that paper for more details. All differential 
equations are integrated in C++ using the adaptive step- 
size integrator StepperDoprS [ii|- First we integrate 
the solution of the homogeneous SN equation ([5]) with 
ingoing boundary conditions at the horizon from = 
2M(1 + Sr) outwards (typically we choose 6r = 10^^). 
Then we integrate the solution with outgoing boundary 



conditions at infinity from Tq 
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typically we choose r, 
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/w inwards, where 
From these two in- 



dependent solutions we compute the Wronskian at the 
(large but finite) radius . We integrate the geodesies 
with given orbital parameters, and at the same time we 
compute the source term corresponding to this trajec- 
tory. We output the solution with ingoing boundary con- 
ditions at the horizon X^^^ and the source term Sim on 
three numerical grids (each consisting of n = 1.6 x 10'^ 
collocation points). The three grids cover the intervals 
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1. Visual comparison of two waveforms for I = m = 2 
with differing cutoffs (e — 10^^, 5 x 10"^). The time axis 
has been shifted so that the maxima of the wave amplitudes 
|^im(i)| (shown in the inset) occur at the same value of {t — 
r*)/M. 



[roo, r+], [nsco, nsco - 0.1] and [nsco - 0.1, r/i], re- 
spectively. This is necessary to keep good accuracy in 
the near-ISCO region, where the particle spends a lot of 
time when e is small. We use a Gauss-Legendre spec- 
tral integrator 11 1 to compute the convolution of X^^'' 



with the source term to find the outgoing wave ampli- 
tude Xi^^{uj). Finally wc sum over multipoles to get 
the total radiated energy, angular momentum and linear 
momentum, and we perform a Fourier transform to get 

The real part of the I = m = 2 waveforms for cutoff 
values e:=5xlO~'^ and e = 10~^ are shown in Figure 
[T] This visual comparison shows quite clearly that the 
ringdown portion of the signal is very weakly dependent 



III. EXTRACTING RINGDOWN AMPLITUDES 
FROM WAVEFORMS 

In this section we describe the procedure we used to ex- 
tract the ringdown amplitudes from the numerical wave- 
forms and to compare with Ref. [l[ . 

In the ringdown phase, the gravitational waveform is 
described by a superposition of quasinormal modes of the 
form 



Xhn = A/" ^ Rnhn exp 



(11) 
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where we dropped the "out" superscript for simplicity. 
Here tOnim arc the (complex) characteristic ringdown fre- 
quencies, Rnim are the amplitudes, Af is an overall nor- 
malization and to is a time shift which determines the 
origin of time. We stress that while some readers may 
be more familiar with the odd and even radiation func- 
tions tplm'^^ (Regge- Wheeler) and V'^m™^ (Zerilh), here 
we work with the Sasaki-Nakamura radiation functions 
Xim, which do not have a specific parity. 
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Given a numerical time-domain waveform Xi,„(i), our 
objective is to extract the amplitudes Rnim- The leading 
amplitude (n = 1) is extracted by plotting log|X;,„| as 
a function of t (see the inset of Figure [IJ . At late times 
only the dominant mode contributes to the signal, and 
log|Xim| becomes a linear function of time with slope 
given by the decay constant 71 = 3(a;ii„i)- The intercept 
is exactly the desired amplitude, log|i?im|- The discrep- 
ancy between the expected and measured 71 contributes 
to the error estimate for ji?!;™]. In practice, this contri- 
bution is minimized if the origin of time is close to the 
center of the linear portion in the curve. 

In the same fashion one can extract the amplitude of 
the second overtone, namely |-R2;m|- In preparation for 
this, one should first extract the phase of Rum- This is 
done by plotting X'^^^ — Xim ■ exp(7iimt) as a function 
of t, so that X' ~ i?exp {i^{ujiim)t) should be a periodic 
function with frequency 3?(a;ii„i) and complex amplitude 
Riim- Once we know the (complex) amplitude Riim, we 
can subtract the dominant ringdown mode to obtain a 
residual X^^{t) = Xi^{t) - Ru ra 0^'piUJilmt. Now we 
can repeat the steps above to find the amplitude of the 
second overtone. Indeed, the procedure can be repeated 
for generic overtone numbers n as long as the signal is 
not dominated by numerical noise, i.e., as long as the 

(n) 

residuals X^^' exhibit a decaying exponential behavior. 

The value of the amplitudes depends, in principle, on 
the constants A/", to- To sidestep this difficulty we com- 
pute the discrepancy of each mode 

A^log(E-^^VKr) , (12) 

where R^^i^^ is the "analytical" amplitude computed in 
Ref. [l| and -R"^™ is the numerically extracted ampli- 
tude. Next we plot A as a function of the mode decay 
constant "fnim- Clearly a rescaling by M will shift A by 
a constant, while a shift of Iq corresponds to a change 
in slope. In particular, for modes with identical decay 
constants (such as modes with different m and fixed I) 
the discrepancy A is the same and it is independent of . 
Now we can carry out a linear fit of A„jm = A„im(7ni)- 
This fit will yield the constants Af, to required for the 
amplitudes to agree, i.e., R^^ = -R""™- The origin of 
time resulting from the linear fit was compared to the 
conventions of Ref. where it what chosen such that 
r{t = 0) = l.lrs. We found that the same happens here 



to a good approximation. 

We performed comparisons for 16 different modes with 
the following {nlm) values (for e = 0.005): the modes 
2 < / < 4 for < m < ? and / = 5 for 2 < m < 5, all 
with n = 1. 

Our main results are shown in Figure [21 After the 
linear fit, the residual errors are seen to be of ~ 1%. We 
also performed the comparison for the mode {nlm) = 
(222) (second overtone of ^ = to = 2). In this case we 
found good agreement up to a residual error of ^ 10%, 
consistent with our error estimate for the numerically 
extracted amplitude in this case. A similar comparison 
was made also for e — 0.01 for some modes, and good 
agreement was found - up to an overall estimated error 
of - 1%. 

In summary, for all the tested modes (see the list 
above) the numerical amplitudes here were found to co- 
incide with the analytic ones computed in [l[. 
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FIG. 2. Amplitude comparison for 16 different modes with 
n = 1 and e = 5 x 10~^. The horizontal axis represents 
the decay constant 7nim in units where rs = 1. Each point 
represents the residual A value (after subtracting a linear fit 
or tuning N ,to) for the given (I, m). Different 7,1;™ values 
correspond to different Vs: the leftmost points correspond to 
I = 2, and 7„i,n grows monotonically with /. 
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